ycls_dat <- read_rds("phc_replications.rds")

fixed <- 1
spacing <- .285

### -------- Study 1: Russian reporters and American news -------------

## Generate pre-covid datasets from published papers
hyman_sheatsley <-
  data.frame(russians_Y = rep(c(0, 1, 0, 1), c(372, 209, 171, 464)),
             russians_Z = rep(c(0, 1), c(581, 635)))

hyman_sheatsley %>%
  group_by(russians_Z) %>%
  summarise(mean(russians_Y))


schuman_presser_1980 <-
  data.frame(russians_Y = rep(c(0, 1, 0, 1), c(155, 187, 85, 250)),
             russians_Z = rep(c(0, 1), c(342, 335)))

schuman_presser_1980 %>%
  group_by(russians_Z) %>%
  summarise(mean(russians_Y))

schuman_presser_1983 <-
  data.frame(russians_Y = rep(c(0, 1, 0, 1), c(65, 52, 32, 75)),
             russians_Z = rep(c(0, 1), c(117, 107)))

schuman_presser_1983 %>%
  group_by(russians_Z) %>%
  summarise(mean(russians_Y))

## Combine with YCLS replication:
combined_df <-
  bind_rows(
    `Hyman and Sheatsley (1950)` = hyman_sheatsley,
    `Schuman and Presser (1980)` = schuman_presser_1980,
    `Schuman et al. (1983)` = schuman_presser_1983,
    `YCLS Week 3` = ycls_dat %>% filter(str_detect(week, "YCLS Week 3")),
    .id = "survey"
  ) %>%
  group_by(survey) %>%
  mutate(russians_Y_s = glass_delta(Y = russians_Y, Z = russians_Z,
                                    level = 0)) 

the_levs <-
  c(
    "Hyman and Sheatsley (1950)",
    "Schuman and Presser (1980)",
    "Schuman et al. (1983)",
    "YCLS Week 3"
  )

the_labs <-
  c(
    "Hyman and Sheatsley (1950)",
    "Schuman and Presser (1980)",
    "Schuman et al. (1983)",
    "Week 3"
  )

## Estimate treatment effects for each sample and combine for plotting
est_study_1 <-
  combined_df %>%
  group_by(survey) %>%
  do(tidy(lm_robust(russians_Y ~ russians_Z, data = .))) %>%
  mutate(estimate_type = "Unstandardized") %>% 
  ## Add standardized estimates for aggregate summary:
  bind_rows(
    combined_df %>%
      group_by(survey) %>%
      do(tidy(lm_robust(russians_Y ~ russians_Z, data = .))) %>% 
      mutate(estimate_type = "Standardized")
  ) %>% 
  filter(term != "(Intercept)") %>%
  mutate(
    survey = factor(
      survey,
      levels = rev(the_levs),
      labels = rev(the_labs)
    ),
    study = ifelse(survey %in% c("Week 3"),
                   "Replication", "Pre-COVID"),
    study = factor(study, levels = c("Replication", "Pre-COVID"))
  )


## Make comparisons
russians_benchmark <- 
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_1 %>% 
      filter(study == "Pre-COVID", 
             estimate_type == "Unstandardized")
  )
russians_benchmark$summary

2*(1-pnorm(russians_benchmark$summary/russians_benchmark$se.summary))

russians_ycls <-
  est_study_1 %>% 
  filter(study == "Replication", 
         estimate_type == "Unstandardized") %>% 
  select(estimate, std.error)

russians_ycls$estimate/russians_benchmark$summary

Z <- (russians_ycls$estimate - russians_benchmark$summary) /
  sqrt((russians_ycls$std.error ^ 2 + russians_benchmark$se.summary ^ 2))
2*(1-pnorm(abs(Z)))

### ----------- Export Figure A.1 for Online Appendix ------
g <-
  est_study_1 %>%
  filter(estimate_type == "Unstandardized") %>%
  ggplot(.,
         aes(
           estimate,
           survey,
           group = study,
           color = study,
           shape = study,
           label = study
         )) +
  geom_vline(xintercept = 0,
             lty = 2,
             size = 0.5) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.5),
    size = 1
  ) +
  geom_point(position = position_dodgev(height = 0.5), size = 2.5) +
  scale_color_grey("", start = 0.2, end = 0.6) +
  scale_shape_manual("", values = c(15, 16)) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  xlab("Average treatment effect estimate") +
  ylab("") +
  ggtitle("") +
  theme_ycls()

g

g %>%
  ggsave(
    filename = "appendix_fig_a1.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * nrow(g$data)
  )


### -------- Study 2: Effect of framing on decision making ---------------

## Re-create original study from published paper
framing_Z = factor(c(rep("Cheap", 93), rep("Expensive", 88)),
                   levels = c("Expensive", "Cheap"))
framing_travel_Y = c(c(rep(1, 63), rep(0, 30)), c(rep(1, 25), rep(0, 63)))
tk_original <-
  data.frame(framing_Z, framing_travel_Y) %>%
  mutate(framing_Z_type = "original")

## Download pre-COVID replications from ManyLabs
tmp <-
  getURL(
    "https://raw.githubusercontent.com/ManyLabsOpenScience/ManyLabs2/master/OSFdata/Framing%20(Tversky%20%26%20Kahneman%2C%201981)/Tversky.1/Global/Data/Tversky_1_study_global_include_all_CLEAN_CASE.csv"
  )

manylabs2_dat <- read.csv(text = tmp) %>%
  mutate(
    framing_travel_Y = as.numeric(coalesce(tver1.1, tver2.1) == 1),
    framing_Z = factor(factor, levels = c("Expensive", "Cheap")),
    framing_Z_type = "original"
  ) %>%
  filter(Setting == "Online (at home)") %>%
  select(framing_travel_Y, framing_Z, framing_Z_type)

## Combine with YCLS replications
combined_df <-
  bind_rows(
    `YCLS Week 7` = ycls_dat %>%
      filter(week == "YCLS Week 7") %>%
      mutate(
        framing_Z  = ifelse(framing_Z == "expensive", "Expensive", "Cheap"),
        framing_Z  = factor(framing_Z, levels = c("Expensive", "Cheap"))
      ) %>%
      select(contains("framing")),
    `Tversky and Kahneman (1981)` = tk_original,
    `Many Labs (2018)` = manylabs2_dat,
    .id = "survey"
  ) %>%
  group_by(survey) %>%
  mutate(framing_travel_Y_s = glass_delta(Y = framing_travel_Y,
                                          Z = framing_Z,
                                          level = "Cheap"))

the_levs <-
  c("Tversky and Kahneman (1981)",
    "Many Labs (2018)",
    "YCLS Week 7")

the_labs <-
  c("Tversky and Kahneman (1981)",
    "Many Labs (2018)",
    "Week 7")

## Estimate treatment effects for each sample and combine for plotting
est_study_2 <-
  combined_df %>%
  group_by(survey, framing_Z_type) %>%
  do(tidy(lm_robust(framing_travel_Y ~ framing_Z, data = .))) %>%
  mutate(estimate_type = "Unstandardized") %>%
  ## Add standardized estimates for aggregate summary:
  bind_rows(
    combined_df %>%
      group_by(survey, framing_Z_type) %>%
      do(tidy(
        lm_robust(framing_travel_Y_s ~ framing_Z, data = .)
      )) %>%
      mutate(estimate_type = "Standardized")
  ) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    survey = factor(
      survey,
      levels = rev(the_levs),
      labels = rev(the_labs)
    ),
    study = ifelse(survey %in% c("Week 7"),
                   "Replication", "Pre-COVID"),
    study = factor(study, levels = c("Replication", "Pre-COVID")),
    type = case_when(
      study == "Replication" &
        framing_Z_type == "modified" ~ "COVID-specific",
      study == "Replication" &
        framing_Z_type == "original" ~ "Direct replication",
      TRUE ~ "Pre-COVID"
    ),
    type = factor(type, levels = rev(
      c("Direct replication", "COVID-specific",
        "Pre-COVID")
    )),
    framing_Z_type = factor(
      framing_Z_type,
      levels = c("original", "modified"),
      labels = c("Original", "Modified")
    )
  )

## Make comparisons
framing_benchmark <- 
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_2 %>% 
      filter(study == "Pre-COVID", 
             estimate_type == "Unstandardized")
  )

framing_benchmark$summary

2*(1-pnorm(framing_benchmark$summary/framing_benchmark$se.summary))

framing_ycls <-
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_2 %>%
      filter(study == "Replication",
             estimate_type == "Unstandardized")
  )

framing_ycls$summary/framing_benchmark$summary

Z <- (framing_ycls$summary - framing_benchmark$summary) /
  sqrt((framing_ycls$se.summary ^ 2 + framing_benchmark$se.summary ^ 2))
2*(1-pnorm(abs(Z)))


### ----------- Export Figure A.2 for Online Appendix ------
g <-
  est_study_2 %>%
  filter(estimate_type == "Unstandardized") %>%
  ggplot(.,
         aes(
           estimate,
           survey,
           color = study,
           shape = study,
           label = type,
           group = type
         )) +
  geom_vline(xintercept = 0,
             lty = 2,
             size = 0.5) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.5),
    size = 1
  ) +
  geom_point(position = position_dodgev(height = 0.5), size = 2.5) +
  scale_color_grey("", start = 0.2, end = 0.6) +
  scale_shape_manual("", values = c(15, 16)) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  xlab("Average treatment effect estimate") +
  ylab("") +
  ggtitle("") +
  theme_ycls() +
  geom_text(
    data = (. %>% filter(survey == "Week 7")),
    aes(x = conf.high + 0.15),
    position = position_dodgev(height = 0.5),
    size = 3,
    family = "serif"
  )

g

g %>%
  ggsave(
    filename = "appendix_fig_a2.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * nrow(g$data)
  )

###-------- Study 3: Gain versus loss framing ---------------

## Many labs replications
manylabs_dat <-
  read_sav("manylabs_asian_disease.sav") %>%
  filter(referrer %in% c("mturk", "pi")) %>%
  mutate(AD_Y_pfp = as.numeric(gainlossDV == 2),
         AD_Z_loss = as.numeric(gainlossgroup == 0)) %>%
  transmute(AD_Y_pfp, AD_Z_loss, referrer)

## Coppock and McClellan + Original
coppock_dat <-
  read_rds("cm_asian_disease.rds") %>%
  mutate(AD_Y_pfp = as.numeric(pfp == 1),
         AD_Z_loss = as.numeric(mortality_frame == 1))

## Combine with YCLS replications
combined_dat <-
  bind_rows(
    `Week 1 - COVID-specific` = ycls_dat %>% filter(week == "YCLS Week 1"),
    `Week 3 - COVID-specific` = ycls_dat %>% filter(week == "YCLS Week 3"),
    `Week 7 - COVID-specific` = ycls_dat %>% filter(week == "YCLS Week 7"),
    `Week 8 - COVID-specific` = ycls_dat %>% filter(week == "YCLS Week 8"),
    `Week 13 - Direct replication` = ycls_dat %>% filter(week == "YCLS Week 13"),
    `Many Labs (2013) - MTurk` = manylabs_dat %>% filter(referrer == "mturk"),
    `Many Labs (2013) - Project Implicit` = manylabs_dat %>% filter(referrer == "pi"),
    `Coppock & McClellan (2016) - MTurk` = coppock_dat %>% filter(survey == "mturk"),
    `Coppock & McClellan (2016) - Lucid` = coppock_dat %>% filter(survey == "lucid"),
    `Kahneman & Tversky (1981)` = coppock_dat %>% filter(survey == "original"),
    .id = "survey"
  ) %>%
  select(starts_with("AD_"), survey) %>%
  ## Exclude arms for Druckman replication in weeks 7, 8, 13
  filter(AD_Z_party_sure_thing == "None" |
           is.na(AD_Z_party_sure_thing),) %>%
  transmute(AD_Y_pfp, AD_Z_loss, survey) %>%
  group_by(survey) %>%
  mutate(AD_Y_pfp_s = glass_delta(Y = AD_Y_pfp, Z = AD_Z_loss,
                                  level = 0)) 

## Estimate treatment effects for each sample and combine for plotting
est_study_3 <-
  combined_dat %>%
  group_by(survey) %>%
  do(tidy(lm_robust(AD_Y_pfp ~ AD_Z_loss, data = .))) %>%
  mutate(estimate_type = "Unstandardized") %>% 
  ## Add standardized estimates for aggregate summary:
  bind_rows(
    combined_dat %>%
      group_by(survey) %>%
      do(tidy(lm_robust(AD_Y_pfp_s ~ AD_Z_loss, data = .))) %>% 
      mutate(estimate_type = "Standardized")
  ) %>% 
  filter(term != "(Intercept)") %>%
  mutate(
    survey = factor(
      survey,
      levels = c(
        "Week 13 - Direct replication",
        "Week 8 - COVID-specific",
        "Week 7 - COVID-specific",
        "Week 3 - COVID-specific",
        "Week 1 - COVID-specific",
        "Coppock & McClellan (2016) - MTurk",
        "Coppock & McClellan (2016) - Lucid",
        "Many Labs (2013) - MTurk",
        "Many Labs (2013) - Project Implicit",
        "Kahneman & Tversky (1981)"
      )
    ),
    study = ifelse(str_detect(survey, "Week"), "Replication", "Pre-COVID"),
    study = factor(study, levels = c("Replication", "Pre-COVID")),
    type = ifelse(survey == "Week 13",
                  "Direct replication",
                  "COVID-specific")
  )

## Make comparisons:
disease_benchmark <- 
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_3 %>% 
      filter(study == "Pre-COVID",
             estimate_type == "Unstandardized")
  )

2*(1-pnorm(disease_benchmark$summary/disease_benchmark$se.summary))

disease_ycls <-
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_3 %>%
      filter(study == "Replication",
             estimate_type == "Unstandardized")
  )

2*(1-pnorm(disease_ycls$summary/disease_ycls$se.summary))

## Compare replications to pre-COVID
disease_ycls$summary/disease_benchmark$summary

Z <- (disease_ycls$summary - disease_benchmark$summary) /
  sqrt((disease_ycls$se.summary ^ 2 + disease_benchmark$se.summary ^ 2))
2*(1-pnorm(abs(Z)))


### ----------- Export Figure A.3 for Online Appendix ------
g <-
  est_study_3 %>%
  filter(estimate_type == "Unstandardized") %>%
  mutate(study_type = paste(study, type)) %>%
  ggplot(.,
         aes(
           estimate,
           survey,
           color = study,
           shape = study_type,
           label = type
         )) +
  geom_vline(xintercept = 0,
             lty = 2,
             size = 0.5) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.5),
    size = 1
  ) +
  geom_point(position = position_dodgev(height = 0.5), size = 2.5) +
  scale_color_grey("", start = 0.2, end = 0.6) +
  scale_shape_manual("", values = c(16, 15, 17)) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  xlab("Average treatment effect estimate") +
  ylab("") +
  ggtitle("") +
  theme_ycls()

g

g %>%
  ggsave(
    filename = "appendix_fig_a3.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * nrow(g$data)
  )

###-------- Study 4: Welfare versus aid to the poor ---------------

## GSS data
gss_dat <- read_rds("gss_welfare.rds") 

## Huber/Paris replications
huber_dat <- read_rds("hp_welfare.rds")

## Combine with YCLS replications.

## First, re-format within-subjects experiments
week_5_atp <-
  ycls_dat %>%
  filter(week == "YCLS Week 5") %>%
  mutate(
    welfare_Y_n = case_when(
      welfare_fedspend_fspr == "Increased" ~ 1,
      welfare_fedspend_fspr == "Kept the same" ~ 0,
      welfare_fedspend_fspr == "Decreased" ~ -1
    ),
    welfare_Z = 0
  ) %>%
  select(welfare_Y_n, welfare_Z)

week_5_welfare <-
  ycls_dat %>%
  filter(week == "YCLS Week 5") %>%
  mutate(
    welfare_Y_n = case_when(
      welfare_fedspend_fswelf == "Increased" ~ 1,
      welfare_fedspend_fswelf == "Kept the same" ~ 0,
      welfare_fedspend_fswelf == "Decreased" ~ -1
    ),
    welfare_Z = 1
  ) %>%
  select(welfare_Y_n, welfare_Z)

week_9_atp <-
  ycls_dat %>%
  filter(week == "YCLS Week 9") %>%
  mutate(
    welfare_Y_n = case_when(
      welfare_fedspend_fspr == "Increased" ~ 1,
      welfare_fedspend_fspr == "Kept the same" ~ 0,
      welfare_fedspend_fspr == "Decreased" ~ -1
    ),
    welfare_Z = 0
  ) %>%
  select(welfare_Y_n, welfare_Z)

week_9_welfare <-
  ycls_dat %>%
  filter(week == "YCLS Week 9") %>%
  mutate(
    welfare_Y_n = case_when(
      welfare_fedspend_fswelf == "Increased" ~ 1,
      welfare_fedspend_fswelf == "Kept the same" ~ 0,
      welfare_fedspend_fswelf == "Decreased" ~ -1
    ),
    welfare_Z = 1
  ) %>%
  select(welfare_Y_n, welfare_Z)

## Treat this as between-subjects experiments by only using response for
## whatever outcome that was randomly assigned to appear first.
week_13 <-
  ycls_dat %>%
  filter(week == "YCLS Week 13") %>%
  mutate(
    welfare_Y_n = case_when(
      welfare_Z_first == 0 ~ welfare_gss_atp_n,
      welfare_Z_first == 1 ~ welfare_gss_welfare_n
    ),
    welfare_Z = welfare_Z_first
  )

combined_dat <-
  bind_rows(
    `Week 1` = ycls_dat %>%
      filter(week == "YCLS Week 1"),
    `Week 2` = ycls_dat %>%
      filter(week == "YCLS Week 2"),
    `Week 3` = ycls_dat %>%
      filter(week == "YCLS Week 3"),
    `Week 4` = ycls_dat %>%
      filter(week == "YCLS Week 4"),
    `Week 5` = bind_rows(week_5_atp, week_5_welfare),
    `Week 6` = ycls_dat %>%
      filter(week == "YCLS Week 6"),
    `Week 7` = ycls_dat %>%
      filter(week == "YCLS Week 7"),
    `Week 8` = ycls_dat %>%
      filter(week == "YCLS Week 8"),
    `Week 9` = bind_rows(week_9_atp, week_9_welfare),
    `Week 11` = ycls_dat %>%
      filter(week == "YCLS Week 11"),
    `Week 12` = ycls_dat %>% 
      filter(week == "YCLS Week 12"),
    `Week 13` = week_13,
    `Huber & Paris (2013) - YouGov` = huber_dat %>%
      filter(survey == "Huber/Paris (MTurk, July '11)"),
    `Huber & Paris (2013) - MTurk` = huber_dat %>%
      filter(survey == "Huber/Paris (YouGov, Apr. '12)"),
    .id = "survey"
  ) %>%
  ## NB: rescale welfare_Z so treatment FX are positive
  bind_rows(gss_dat) %>% 
  group_by(survey) %>%
  mutate(
    welfare_Z = as.numeric(welfare_Z == 0),
    welfare_Y_s = glass_delta(Y = welfare_Y_n, Z = welfare_Z,
                              level = 0)
  ) %>%
  select(welfare_Y_n, welfare_Y_s, welfare_Z, survey) 

## Label for ordering of estimates:
gss_levs <-
  combined_dat %>%
  filter(str_detect(survey, 'GSS')) %>%
  pull(survey) %>%
  unique()

lucid_levs <-
  combined_dat %>%
  filter(str_detect(survey, 'Week')) %>%
  pull(survey) %>%
  unique()

the_levs <-
  c(rev(gss_levs),
    rev(lucid_levs),
    "Huber & Paris (2013) - YouGov",
    "Huber & Paris (2013) - MTurk")

## Estimate treatment effects for each sample and combine for plotting
est_study_4 <-
  combined_dat %>% 
  group_by(survey) %>%
  do(tidy(lm_robust(welfare_Y_s ~ welfare_Z, data = .))) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    survey = factor(survey, levels = the_levs),
    mode = case_when(
      survey %in% paste("GSS", 2002:2018)  ~ "Computer-assisted personal interviews",
      survey %in% paste("GSS", 1986:2000) ~ "Paper and pencil interviews",
      TRUE ~ "Computer-assisted web interviews"
    ),
    mode = factor(
      mode,
      levels = c(
        "Computer-assisted web interviews",
        "Computer-assisted personal interviews",
        "Paper and pencil interviews"
      )
    ),
    study = case_when(
      mode %in% c(
        "Computer-assisted web interviews",
        "Paper and pencil interviews"
      ) &
        str_detect(survey, "Huber") ~ "Pre-COVID",
      mode == "Computer-assisted web interviews" &
        !str_detect(survey, "Huber") ~ "Replication",
      mode != "Computer-assisted web interviews" ~ "Pre-COVID"
    ),
    study = factor(study, levels = c("Replication", "Pre-COVID")),
    type = ifelse(
      survey %in% c("Week 5", "Week 9"),
      "Observational",
      "Experimental"
    ),
    estimate_type = "Standardized"
  )

## Make comparisons 
ycls_cawi <- 
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_4 %>% 
      filter(mode == "Computer-assisted web interviews",
             type == "Experimental", 
             estimate_type == "Standardized")
  )

2*(1-pnorm(ycls_cawi$summary/ycls_cawi$se.summary))

## GSS CAPI estimates summary
gss_capi <-
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_4 %>%
      filter(
        str_detect(survey, "GSS"),
        type == "Experimental",
        mode == "Computer-assisted personal interviews",
        estimate_type == "Standardized"
      )
  )

gss_capi$summary
gss_capi$se.summary

## GSS PAPI estimates summary
gss_papi <-
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_4 %>%
      filter(
        str_detect(survey, "GSS"),
        type == "Experimental",
        mode == "Paper and pencil interviews",
        estimate_type == "Standardized"
      )
  )

gss_papi$summary
gss_papi$se.summary


## Experimental YCLS estimates summary
ycls_exp <- 
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_4 %>% 
      filter(study == "Replication", 
             type == "Experimental",
             estimate_type == "Standardized")
  )

2*(1-pnorm(ycls_exp$summary/ycls_exp$se.summary))

ycls_exp$se.summary
ycls_exp$summary

## Pre-COVID CAWI benchmark
huber_benchmark <-
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_4 %>%
      filter(str_detect(survey, "Huber"),
             estimate_type == "Standardized")
  )

huber_benchmark$se.summary
huber_benchmark$summary

## Pre-COVID benchmarks (Huber + Original GSS)
ycls_benchmark <- 
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_4  %>% 
      filter(str_detect(survey, "Huber") | survey == "GSS 1986",
             estimate_type == "Standardized")
  )

ycls_benchmark$se.summary
ycls_benchmark$summary

## Observational YCLS estimates summary
ycls_obs <-
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_4 %>%
      filter(
        study == "Replication",
        type != "Experimental",
        estimate_type == "Standardized"
      )
  )

ycls_obs$se.summary
ycls_obs$summary

### ----------- Export Figure A.4 for Online Appendix ------
g <-
  est_study_4 %>%
  mutate(study_type = paste(study, type)) %>%
  ggplot(.,
         aes(
           estimate,
           survey,
           group = type,
           color = study,
           shape = study_type,
           label = type
         )) +
  geom_vline(xintercept = 0,
             lty = 2,
             size = 0.5) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.5),
    size = 1
  ) +
  geom_point(position = position_dodgev(height = 0.5), size = 2.5) +
  facet_col(~ mode, scales = "free_y", space = "free") +
  scale_x_continuous(breaks = round(seq(-1.4, 1.2, by = .2), 2)) +
  scale_color_grey("", start = 0.2, end = 0.6) +
  scale_shape_manual("", values = c(16, 15, 17)) +
  xlab("Average treatment effect estimate (standard units)") +
  ylab("") +
  ggtitle("") +
  theme_ycls() +
  geom_text(
    data = (. %>% filter(survey %in% c("Week 9", "Week 11"))),
    aes(x = conf.high + 0.2),
    position = position_dodgev(height = 0.5),
    size = 3,
    family = "serif"
  )

g

g %>%
  ggsave(
    filename = "appendix_fig_a4.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * (nrow(g$data) * 0.75)
  )

###-------- Study 5: Gain versus loss framing + party endorsements ------------

## Data from original study
druckman <- 
  read_csv("druckman_original.csv") %>% 
  rename(x_pid3 = partyid)

## Combine with YCLS replications
combined_dat <-
  bind_rows(
    `Week 7` = ycls_dat %>% 
      filter(week == "YCLS Week 7"),
    `Week 8` = ycls_dat %>% 
      filter(week == "YCLS Week 8"),
    `Week 13 (direct replication)` = ycls_dat %>% 
      filter(week == "YCLS Week 13"),
    `Druckman (2001)` = druckman,
    .id = "survey"
  ) %>%
  select(starts_with("AD_"), x_pid3, survey) %>% 
  mutate(
    AD_Z_party_sure_thing = factor(
      AD_Z_party_sure_thing,
      levels = rev(c("None", "Democrats", "Republicans")),
      labels = rev(c("Program A", "Democrats'\n Program", "Republicans'\n Program"))
    ),
    x_pid3 = factor(
      x_pid3,
      levels = c("Democrat", "Independent", "Republican"),
      labels = c("Framing effects by label of risk-averse alternative, among Democrats", 
                 "Framing effects by label of risk-averse alternative, among Independents", 
                 "Framing effects by label of risk-averse alternative, among Republicans")
    )
  ) %>%
  filter(!is.na(x_pid3),!is.na(AD_Z_party_sure_thing))

## Estimate treatment effects for each sample and combine for plotting
est_study_5 <-
  combined_dat %>%
  group_by(survey, x_pid3, AD_Z_party_sure_thing) %>%
  do(tidy(lm_robust(AD_Y_sure_thing ~ AD_Z_gain, data = .))) %>%
  mutate(estimate_type = "Unstandardized") %>%
  ## Add standardized estimates for aggregate summary:
  bind_rows(
    combined_dat %>%
      group_by(survey, AD_Z_party_sure_thing) %>%
      mutate(AD_Y_sure_thing_s = glass_delta(
        Y = AD_Y_sure_thing, Z = AD_Z_gain,
        level = 0
      )) %>%
      do(tidy(
        lm_robust(AD_Y_sure_thing_s ~ AD_Z_gain, data = .)
      )) %>%
      mutate(estimate_type = "Standardized")
  ) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    study = ifelse(str_detect(survey, "Week"), "Replication", "Pre-COVID"),
    study = factor(study, levels = c("Replication", "Pre-COVID")),
    type = ifelse(
      survey == "Week 13 (direct replication)",
      "Direct replication",
      "COVID-specific"
    ),
    survey = factor(survey, levels = rev(
      c(
        "Druckman (2001)",
        "Week 7",
        "Week 8",
        "Week 13 (direct replication)"
      )
    ))
  ) 

## Make comparisons
druckman_benchmark <-
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_5 %>%
      filter(study == "Pre-COVID",
             estimate_type == "Unstandardized")
  )

2*(1-pnorm(druckman_benchmark$summary/druckman_benchmark$se.summary))

druckman_ycls <- 
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_5 %>%
      filter(study == "Replication",
             estimate_type == "Unstandardized")
  )

2*(1-pnorm(druckman_ycls$summary/druckman_ycls$se.summary))

## Compare replications to pre-COVID
druckman_ycls$summary/druckman_benchmark$summary

Z <- (druckman_ycls$summary - druckman_benchmark$summary) /
  sqrt((druckman_ycls$se.summary ^ 2 + druckman_benchmark$se.summary ^ 2))
2*(1-pnorm(abs(Z)))


### ----------- Export Figure A.5 for Online Appendix ------
g <-
  est_study_5 %>%
  filter(x_pid3 != "Framing effects by label of risk-averse alternative, among Independents",
         estimate_type != "Standardized") %>%
  ggplot(.,
         aes(
           estimate,
           AD_Z_party_sure_thing,
           group = survey,
           color = study,
           shape = study,
           label = survey
         )) +
  geom_vline(xintercept = 0,
             lty = 2,
             size = 0.5) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(0.5),
    size = 1
  ) +
  geom_point(position = position_dodgev(0.5),
             size = 2.5,
             fill = "black") +
  scale_shape_manual("", values = c(15, 16)) +
  scale_color_grey("", start = 0.2, end = 0.6) +
  facet_wrap(~ x_pid3, ncol = 1) +
  coord_cartesian(xlim = c(-0.45, 0.85)) +
  scale_x_continuous(labels = scales::percent,
                     breaks = seq(-0.4, 0.8, by = 0.2)) +
  xlab("Average treatment effect estimate") +
  ylab("") +
  ggtitle("") +
  theme_ycls() +
  geom_text(
    data = (
      . %>% filter(
        AD_Z_party_sure_thing == "Program A",
        x_pid3 == "Framing effects by label of risk-averse alternative, among Democrats"
      )
    ),
    aes(label = survey, x = c(.80, .60, .35, .45)),
    size = 2.5,
    position = position_dodgev(0.5)
  )


g


g %>%
  ggsave(
    filename = "appendix_fig_a5.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * (nrow(g$data)*.90)
  )

### ----------- Export Table A.5 for Online Appendix ------

## Original results x endorsement 
druckman_estimates <-
  est_study_5 %>%
  filter(
    str_detect(survey, "Druckman"),
    !str_detect(x_pid3, "Independent"),
    estimate_type == "Unstandardized"
  ) %>%
  mutate(
    AD_Z_party_sure_thing = case_when(
      str_detect(AD_Z_party_sure_thing, "Program A") ~ "Program A",
      str_detect(AD_Z_party_sure_thing, "Republican") ~ "Republicans' Program",
      str_detect(AD_Z_party_sure_thing, "Democrat") ~ "Democrats' Program"
    ),
    x_pid3 = case_when(
      str_detect(x_pid3, "Republican") ~ "Republicans",
      str_detect(x_pid3, "Democrat") ~ "Democrats"
    )
  ) %>%
  select(survey, AD_Z_party_sure_thing, x_pid3, estimate, std.error) 

ycls_estimates <-
  est_study_5 %>%
  filter(
    !str_detect(survey, "Druckman"),!str_detect(x_pid3, "Independent"),
    estimate_type == "Unstandardized"
  ) %>%
  group_by(x_pid3, AD_Z_party_sure_thing) %>%
  summarise(estimate = weighted.mean(estimate, w = 1 / (std.error) ^ 2),
            std.error = sqrt(1 / sum(1 / (std.error) ^ 2))) %>%
  mutate(
    AD_Z_party_sure_thing = case_when(
      str_detect(AD_Z_party_sure_thing, "Program A") ~ "Program A",
      str_detect(AD_Z_party_sure_thing, "Republican") ~ "Republicans' Program",
      str_detect(AD_Z_party_sure_thing, "Democrat") ~ "Democrats' Program"
    ),
    x_pid3 = case_when(
      str_detect(x_pid3, "Republican") ~ "Republicans",
      str_detect(x_pid3, "Democrat") ~ "Democrats"
    )
  ) %>%
  select(AD_Z_party_sure_thing, x_pid3, estimate, std.error) %>%
  mutate(survey = "YCLS summary") %>%
  select(survey, everything())

combined_estimates <- 
  bind_rows(druckman_estimates, ycls_estimates) %>% 
  arrange(x_pid3, AD_Z_party_sure_thing) %>% 
  select(AD_Z_party_sure_thing, estimate, x_pid3) %>% 
  spread(survey, estimate) %>% 
  mutate(Difference = `YCLS summary` - `Druckman (2001)`,
         Ratio = `YCLS summary`/`Druckman (2001)`) %>% 
  gather(survey, estimate, `Druckman (2001)`:Difference) 

combined_ses <- 
  bind_rows(druckman_estimates, ycls_estimates) %>% 
  arrange(x_pid3, AD_Z_party_sure_thing) %>% 
  select(AD_Z_party_sure_thing, std.error, x_pid3) %>% 
  spread(survey, std.error) %>% 
  mutate(Difference = sqrt(`Druckman (2001)`^2 + `YCLS summary`^2)) %>% 
  gather(survey, std.error, `Druckman (2001)`:Difference) 

## Table of output 
druckman_summary <- 
  left_join(combined_estimates, combined_ses) %>% 
  mutate(p = 2*(1-pnorm(abs(estimate/std.error))), 
         entry = make_entry(est = estimate, se = std.error, p = p, digits = 2),
         AD_Z_party_sure_thing = factor(AD_Z_party_sure_thing, levels = c("Program A", "Democrats' Program",
                                          "Republicans' Program"))) %>% 
  select(AD_Z_party_sure_thing, entry, Ratio, x_pid3, survey) %>%
  spread(survey, entry) %>% 
  select(AD_Z_party_sure_thing, x_pid3, `YCLS summary`,  `Druckman (2001)`, 
         Difference, Ratio) 

xtable(druckman_summary) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0("appendix_table_a5.tex"))


###-------- Study 6: Foreign aid misperceptions ------------
gilens <- read_rds("gilens_original.rds")

combined_dat <-
  bind_rows(`Gilens (2001)` = gilens,
            `YCLS Week 3` = ycls_dat %>% 
              filter(week == "YCLS Week 3"),
            .id = "survey") %>% 
  group_by(survey) %>% 
  mutate(gilens_Y_std = glass_delta(Y = gilens_Y, Z = gilens_Z,
                                    level = 0))

the_levs <-
  c("Gilens (2001)",
    "YCLS Week 3")

the_labs <- 
  c("Gilens (2001)",
    "Week 3")

## Estimate treatment effects for each sample and combine for plotting
est_study_6 <-
  combined_dat %>%
  group_by(survey) %>%
  do(tidy(lm_robust(gilens_Y_std ~ gilens_Z, data = .))) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    survey = factor(
      survey,
      levels = rev(the_levs),
      labels = rev(the_labs)
    ),
    study = ifelse(survey %in% c("Week 3"),
                   "Replication", "Pre-COVID"),
    study = factor(study, levels = c("Replication", "Pre-COVID")),
    estimate_type = "Standardized"
  )

## Make comparisons
ycls_orig <-
  est_study_6 %>% 
  filter(survey != "Week 3")

ycls_replication <- 
  est_study_6 %>% 
  filter(survey == "Week 3")

Z <- (ycls_orig$estimate - ycls_replication$estimate)/
  (sqrt(ycls_orig$std.error^2 + ycls_replication$std.error^2))
2*(1-pnorm(abs(Z)))


### ----------- Export Figure A.6 for Online Appendix ------
g <-
  ggplot(est_study_6, aes(estimate, survey, group = study, 
                    color = study, shape = study, 
                    label = study)) +
  geom_vline(xintercept = 0,
             lty = 2,
             size = 0.5) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.5),
    size = 1
  ) +
  geom_point(position = position_dodgev(height = 0.5), size = 2.5) +
  scale_color_grey("", start = 0.2, end = 0.6) +
  scale_shape_manual("", values = c(15, 16)) +
  xlab("Average treatment effect estimate (standard units)") +
  ylab("") +
  ggtitle("") + 
  theme_ycls() 

g

g %>%
  ggsave(
    filename = "appendix_fig_a6.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * nrow(g$data)
  )

###-------- Study 7: Perceived intentionality for side effects ------------

## Download pre-COVID replications from ManyLabs
tmp <-
  getURL(
    "https://raw.githubusercontent.com/ManyLabsOpenScience/ManyLabs2/master/OSFdata/Intentional%20Side-Effects%20(Knobe%2C%202003)/Knobe.1/Global/Data/Knobe_1_study_global_include_all_CLEAN_CASE.csv"
  )

manylabs2_dat <- read.csv(text = tmp) %>%
  mutate(
    knobe_intent_Y_n = coalesce(knob1.3, knob2.3),
    knobe_Z = factor(factor, levels = c("Help", "Harm")),
    knobe_intent_binary = as.numeric(knobe_intent_Y_n > 4),
    knobe_Z_type = "original"
  ) %>%
  filter(Setting == "Online (at home)")


## Re-create original study
knobe_Z = factor(c(rep("Help", 39), rep("Harm", 39)),
                 levels = c("Help", "Harm"))
knobe_intent_binary = c(c(rep(1, 9), rep(0, 30)), c(rep(1, 32), rep(0, 7)))
knobe_original <-
  data.frame(knobe_Z, knobe_intent_binary) %>%
  mutate(knobe_Z_type = "original")

## Week 7 replication
replication <- 
  ycls_dat %>% 
  filter(week == "YCLS Week 7") %>% 
  mutate(knobe_Z = ifelse(knobe_Z == "harmed", "Harm", "Help"),
         knobe_Z = factor(knobe_Z, levels = c("Help", "Harm")),
         knobe_intent_binary = as.numeric(knobe_intent_Y_n > 4)) %>% 
  select(starts_with("knobe_"))

## Combine datsets:
combined_dat <-
  bind_rows(
    `Week 7` = replication,
    `Knobe (2003)` = knobe_original,
    `Many Labs (2018)` = manylabs2_dat,
    .id = "survey"
  ) %>%
  group_by(survey, knobe_Z_type) %>%
  mutate(knobe_intent_binary_s = glass_delta(Y = knobe_intent_binary,
                                             Z = knobe_Z, level = "Help"))

the_labs <-
  c("Knobe (2003)",
    "Many Labs (2018)",
    "Week 7")

## Estimate treatment effects for each sample and combine for plotting
est_study_7 <-
  combined_dat %>%
  group_by(survey, knobe_Z_type) %>%
  do(tidy(lm_robust(knobe_intent_binary ~ knobe_Z, data = .))) %>%
  mutate(estimate_type = "Unstandardized") %>% 
  ## Add standardized estimates for aggregate summary:
  bind_rows(
    combined_dat %>%
      group_by(survey, knobe_Z_type) %>%
      do(tidy(lm_robust(knobe_intent_binary_s ~ knobe_Z, data = .))) %>%
      mutate(estimate_type = "Standardized") 
  ) %>% 
  filter(term != "(Intercept)") %>%
  mutate(
    survey = factor(survey, levels = rev(the_labs)),
    study = ifelse(survey %in% c("Week 7"), 
                   "Replication", "Pre-COVID"),
    study = factor(study, levels = c("Replication", "Pre-COVID")),
    type = case_when(
      study == "Replication" & knobe_Z_type == "modified" ~ "COVID-specific",
      study == "Replication" & knobe_Z_type == "original" ~ "Direct replication",
      TRUE ~ "Pre-COVID"
    ),
    type = factor(type, levels = c("Direct replication", "COVID-specific", 
                                   "Pre-COVID")), 
    knobe_Z_type = factor(
      knobe_Z_type,
      levels = c("original", "modified"),
      labels = c("Original", "Modified")
    )
  )


## Comparisons
knobe_benchmark <- 
  rmeta::meta.summaries(
    d = estimate,
    se = std.error,
    data = est_study_7 %>% 
      filter(survey != "Week 7", estimate_type == "Unstandardized")
  )

ycls_direct <- 
  est_study_7 %>% 
  filter(survey == "Week 7", type == "Direct replication",
         estimate_type == "Unstandardized")

ycls_direct$estimate/knobe_benchmark$summary

Z <- (knobe_benchmark$summary - ycls_direct$estimate)/
  (sqrt(knobe_benchmark$se.summary^2 + ycls_direct$std.error^2))

2*(1-pnorm(abs(Z)))

## Compare YCLS direct w/ covid-specific
ycls_covid <- 
  est_study_7 %>% 
  filter(survey == "Week 7", type != "Direct replication",
         estimate_type == "Unstandardized")

Z <- (ycls_covid$estimate - ycls_direct$estimate) /
  (sqrt(ycls_covid$std.error ^ 2 + ycls_direct$std.error ^ 2))

2*(1-pnorm(abs(Z)))

### ----------- Export Figure A.7 for Online Appendix ------
g <-
  est_study_7 %>% 
  filter(estimate_type == "Unstandardized") %>% 
  mutate(
    survey = case_when(
      survey == "Week 7" & knobe_Z_type == "Original" ~ "Direct replication",
      survey == "Week 7" & knobe_Z_type == "Modified" ~ "COVID-specific",
      TRUE ~ paste(survey)
    ),
    survey = factor(survey, levels = rev(c("Knobe (2003)",
                                           "Many Labs (2018)",
                                           "Direct replication",
                                           "COVID-specific")))
  ) %>% 
  ggplot(., aes(estimate, survey, 
                color = study, 
                shape = type)) +
  geom_vline(xintercept = 0,
             lty = 2,
             size = 0.5) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.5),
    size = 1
  ) +
  geom_point(position = position_dodgev(height = 0.5), size = 2.5) +
  scale_color_grey("", start = 0.2, end = 0.6) +
  scale_shape_manual("", values = c(15, 17, 16)) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  xlab("Average treatment effect estimate") +
  ylab("") +
  ggtitle("") + 
  theme_ycls()

g

g %>%
  ggsave(
    filename = "appendix_fig_a7.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * nrow(g$data)
  )

###-------- Study 8: Atomic aversion ------------

## Combine original and Aronow et al. replication w/ YCLS estimates
psv_dat <- read_rds("psv_original.rds")
abp_dat <- read_rds("abp_replication.rds")

combined_dat <-
  bind_rows(
    `Press et al. (2013)` = psv_dat,
    `Aronow et al. (2019)` = abp_dat,
    `Week 5` = ycls_dat %>%
      filter(week == "YCLS Week 5"),
    `Week 6` = ycls_dat %>%
      filter(week == "YCLS Week 6"),
    `Week 13` = ycls_dat %>%
      filter(week == "YCLS Week 13"),
    .id = "survey"
  ) %>%
  filter(is.na(psv_Z_video) | psv_Z_video == 0) %>%
  mutate(
    survey = factor(
      survey,
      levels = c(
        "Week 13",
        "Week 6",
        "Week 5",
        "Aronow et al. (2019)",
        "Press et al. (2013)"
      )
    ),
    psv_Z_prospective = factor(
      psv_Z_prospective,
      levels = c(
        "Nukes 90%; Conventional 90%",
        "Nukes 90%; Conventional 70%",
        "Nukes 90%; Conventional 45%"
      )
    ),
    psv_approve_nuke_binary_s = glass_delta(Y = psv_approve_nuke_binary,
                                            Z = psv_Z_prospective,
                                            level = "Nukes 90%; Conventional 90%"),
    psv_prefer_binary_s = glass_delta(Y = psv_prefer_binary,
                                      Z = psv_Z_prospective,
                                      level = "Nukes 90%; Conventional 90%"),
    psv_retrospective_approve_binary_s = glass_delta(Y = psv_retrospective_approve_binary,
                                                     Z = psv_Z_retrospective,
                                                     level = "Conventional Strike"),
    psv_retrospective_ethical_binary_s = glass_delta(Y = psv_retrospective_ethical_binary,
                                                     Z = psv_Z_retrospective,
                                                     level = "Conventional Strike")
    
  ) 
  

the_levs <-
  c("Week 13",
    "Week 6",
    "Week 5",
    "Aronow et al. (2019)",
    "Press et al. (2013)")

## Estimate treatment effects for each sample and combine for plotting

## Prospective treatment/outcome pairs
approve_effects <-
  combined_dat %>%
  group_by(survey) %>%
  do(tidy(
    lm_robust(psv_approve_nuke_binary ~ psv_Z_prospective, data = .)
  )) %>% 
  mutate(estimate_type = "Unstandardized") %>% 
  ## Add standardized estimates for aggregate summary:
  bind_rows(
    combined_dat %>%
      group_by(survey) %>%
      do(tidy(
        lm_robust(psv_approve_nuke_binary_s ~ psv_Z_prospective, data = .)
      )) %>% 
      mutate(estimate_type = "Standardized") 
  ) 

prefer_effects <-
  combined_dat %>%
  group_by(survey) %>%
  do(tidy(lm_robust(
    psv_prefer_binary  ~ psv_Z_prospective, data = .
  ))) %>% 
  mutate(estimate_type = "Unstandardized") %>% 
  ## Add standardized estimates for aggregate summary:
  bind_rows(
    combined_dat %>%
      group_by(survey) %>%
      do(tidy(
        lm_robust(psv_prefer_binary_s ~ psv_Z_prospective, data = .)
      )) %>% 
      mutate(estimate_type = "Standardized") 
  ) 


est_study_8a <-
  bind_rows(approve_effects,
            prefer_effects) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    term = gsub(pattern = "psv_Z_prospective", "", term),
    outcome_group =
      case_when(
        str_detect(outcome, "psv_approve") ~ "Approve Nuclear Use",
        str_detect(outcome, "psv_prefer")  ~ "Prefer Nuclear Use"
      ),
    outcome_group = factor(
      outcome_group,
      levels = c("Prefer Nuclear Use",
                 "Approve Nuclear Use")
    ),
    survey = factor(survey, levels = the_levs),
    study = ifelse(
      survey %in% c("Week 5", "Week 6", "Week 13"),
      "Replication",
      "Pre-COVID"
    ),
    study = factor(study, levels = c("Replication", "Pre-COVID")),
    term = factor(
      term,
      levels = c("Nukes 90%; Conventional 70%",
                 "Nukes 90%; Conventional 45%"),
      labels = c("90/70",
                 "90/45")
    ),
    term_group = case_when(
      term == "Nukes 90%; Conventional 70%" &
        study == "Pre-COVID" ~ "type1",
      term == "Nukes 90%; Conventional 70%" &
        study == "Replication" ~ "type2",
      term == "Nukes 90%; Conventional 45%" ~ "type3"
    )
  )

## Retrospective treatment/outcome pairs
approve_effects <-
  combined_dat %>%
  group_by(survey) %>%
  do(tidy(
    lm_robust(psv_retrospective_approve_binary ~ psv_Z_retrospective, data = .)
  )) %>% 
  mutate(estimate_type = "Unstandardized") %>% 
  ## Add standardized estimates for aggregate summary:
  bind_rows(
    combined_dat %>%
      group_by(survey) %>%
      do(tidy(
        lm_robust(psv_retrospective_approve_binary_s ~ psv_Z_retrospective, 
                  data = .)
      )) %>% 
      mutate(estimate_type = "Standardized") 
  ) 


ethical_effects <-
  combined_dat %>%
  group_by(survey) %>%
  do(tidy(
    lm_robust(psv_retrospective_ethical_binary  ~ psv_Z_retrospective, data = .)
  )) %>% 
  mutate(estimate_type = "Unstandardized") %>% 
  ## Add standardized estimates for aggregate summary:
  bind_rows(
    combined_dat %>%
      group_by(survey) %>%
      do(tidy(
        lm_robust(psv_retrospective_ethical_binary_s  ~ psv_Z_retrospective, 
                  data = .)
      )) %>% 
      mutate(estimate_type = "Standardized") 
  ) 

est_study_8b <-
  bind_rows(approve_effects,
            ethical_effects) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    term = "Nuclear",
    outcome_group =
      case_when(
        outcome == "psv_retrospective_approve_binary" ~ "Approve Strike",
        outcome == "psv_retrospective_ethical_binary" ~ "Ethical Strike"
      ),
    survey = factor(survey, levels = the_levs),
    study = ifelse(
      survey %in% c("Week 5", "Week 6", "Week 13"),
      "Replication",
      "Pre-COVID"
    ),
    study = factor(study, levels = c("Replication", "Pre-COVID"))
  )


### ----------- Export Figure A.8a for Online Appendix ------
g <- 
  est_study_8a %>% 
  filter(estimate_type == "Unstandardized") %>% 
  ggplot(.,
         aes(
           estimate,
           survey,
           group = term,
           color = study,
           shape = term
         )) +
  geom_vline(
    xintercept = 0,
    lty = 2,
    size = 0.5
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.5), size = 1
  ) +
  geom_point(position = position_dodgev(height = 0.5), size = 2.5 
  ) +
  facet_wrap(~ outcome_group) + 
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_color_grey("", start = 0.2, end = 0.6) +
  xlab("Average treatment effect estimate") +
  ylab("") +
  theme_ycls() + 
  theme(axis.line.x = element_line(size = 0.5),
        strip.text.x = element_text(hjust = 0.5)) + 
  coord_capped_cart(bottom = "none") +
  geom_text(
    data = (. %>% filter(survey ==  "Press et al. (2013)", 
                         outcome_group == "Prefer Nuclear Use")),
    aes(label = term, x = conf.low - 0.1), color = "black",
    position = position_dodgev(height = 0.5),
    size = 3, family = "serif") 
g


g %>%
  ggsave(
    filename = "appendix_fig_a8a.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * (nrow(g$data)*0.8)
  )

### ----------- Export Table A.8a for Online Appendix ------

psv_prospective_estimates <- 
  est_study_8a %>% 
  ungroup() %>% 
  filter(survey == "Press et al. (2013)",
         estimate_type == "Unstandardized") %>% 
  mutate(Z_psv = term,
         psv_type = "Prospective") %>%
  select(survey, Z_psv, psv_type, outcome_group, estimate, std.error) 

abp_prospective_estimates <- 
  est_study_8a %>% 
  ungroup() %>% 
  filter(survey == "Aronow et al. (2019)",
         estimate_type == "Unstandardized") %>% 
  mutate(Z_psv = term,
         psv_type = "Prospective") %>%
  select(survey, Z_psv, psv_type, outcome_group, estimate, std.error) 

ycls_prospective_estimates <- 
  est_study_8a %>% 
  ungroup() %>% 
  filter(str_detect(survey, "Week")) %>% 
  mutate(Z_psv = term) %>%
  group_by(Z_psv, outcome_group) %>% 
  summarise(estimate = weighted.mean(estimate, w = 1/(std.error)^2),
            std.error = sqrt(1 / sum(1 / (std.error) ^ 2))) %>% 
  mutate(survey = "YCLS summary",
         psv_type = "Prospective") %>% 
  select(survey, everything())

combined_prospective_estimates <- 
  bind_rows(psv_prospective_estimates, abp_prospective_estimates, 
            ycls_prospective_estimates) %>% 
  select(survey, Z_psv, outcome_group, estimate) %>% 
  pivot_wider(names_from = survey, values_from = estimate) %>% 
  mutate(abp_diff = `YCLS summary` - `Aronow et al. (2019)`,
         psv_diff = `YCLS summary` - `Press et al. (2013)`,
         abp_ratio = `YCLS summary`/`Aronow et al. (2019)`,
         psv_ratio = `YCLS summary`/`Press et al. (2013)`) %>% 
  gather(survey, estimate, `Press et al. (2013)`:psv_diff) 

combined_prospective_ses <- 
  bind_rows(psv_prospective_estimates, abp_prospective_estimates, 
            ycls_prospective_estimates) %>% 
  select(survey, Z_psv, outcome_group, std.error) %>% 
  pivot_wider(names_from = survey, values_from = std.error) %>% 
  mutate(abp_diff = sqrt(`Aronow et al. (2019)`^2 + `YCLS summary`^2),
         psv_diff = sqrt(`Press et al. (2013)`^2 + `YCLS summary`^2)) %>% 
  gather(survey, std.error, `Press et al. (2013)`:psv_diff) 

## Table of output 
prospective_summary <- 
  left_join(combined_prospective_estimates, combined_prospective_ses) %>% 
  mutate(p = 2*(1-pnorm(abs(estimate/std.error))), 
         entry = make_entry(est = estimate, se = std.error, p = p, digits = 2),
         outcome_group = ifelse(outcome_group == "Prefer Nuclear Use",
                                "Prefer", "Approve"),
         abp_ratio = ifelse(abp_ratio < 0, "-", 
                            paste0(sprintf(paste0("%.", 2, "f"), abp_ratio))),
         psv_ratio = ifelse(psv_ratio < 0, "-", 
                            paste0(sprintf(paste0("%.", 2, "f"), psv_ratio)))) %>% 
  select(Z_psv, entry, abp_ratio, psv_ratio, outcome_group, survey) %>%
  spread(survey, entry) %>% 
  select(Z_psv, outcome_group, `YCLS summary`, `Press et al. (2013)`, psv_diff, 
         psv_ratio, `Aronow et al. (2019)`, abp_diff, abp_ratio) %>% 
  arrange(outcome_group) 

prospective_summary %>% 
  filter(outcome_group == "Prefer") %>% 
  xtable() %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0("appendix_table_a8a_top.tex"))

prospective_summary %>% 
  filter(outcome_group == "Approve") %>% 
  xtable() %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0("appendix_table_a8a_btm.tex"))

### ----------- Export Figure A.8b for Online Appendix ------
g <- 
  est_study_8b %>% 
  filter(estimate_type == "Unstandardized") %>% 
  ggplot(.,
         aes(
           estimate,
           survey,
           color = study,
           shape = term
         )) +
  geom_vline(
    xintercept = 0,
    lty = 2,
    size = 0.5
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.5), size = 1
  ) +
  geom_point(position = position_dodgev(height = 0.5), size = 2 
  ) +
  facet_wrap(~ outcome_group) + 
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_color_grey("", start = 0.2, end = 0.6) +
  xlab("Average treatment effect estimate") +
  ylab("") +
  theme_ycls() + 
  theme(axis.line.x = element_line(size = 0.5),
        strip.text.x = element_text(hjust = 0.5)) + 
  coord_capped_cart(bottom = "none") 

g

g %>%
  ggsave(
    filename = "appendix_fig_a8b.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * (nrow(g$data))
  )

### ----------- Export Table A.8b for Online Appendix ------
psv_retrospective_estimates <- 
  est_study_8b %>% 
  ungroup() %>% 
  filter(survey == "Press et al. (2013)",
         estimate_type == "Unstandardized") %>% 
  mutate(Z_psv = term, 
         psv_type = "retrospective") %>%
  select(survey, Z_psv, psv_type, outcome_group, estimate, std.error) 

abp_retrospective_estimates <- 
  est_study_8b %>% 
  ungroup() %>% 
  filter(survey == "Aronow et al. (2019)",
         estimate_type == "Unstandardized") %>% 
  mutate(Z_psv = term, 
         psv_type = "retrospective") %>%
  select(survey, Z_psv, psv_type, outcome_group, estimate, std.error) 

ycls_retrospective_estimates <- 
  est_study_8b %>% 
  ungroup() %>% 
  filter(str_detect(survey, "Week"),
         estimate_type == "Unstandardized") %>% 
  mutate(Z_psv = term) %>%
  group_by(Z_psv, outcome_group) %>% 
  summarise(estimate = weighted.mean(estimate, w = 1/(std.error)^2),
            std.error = sqrt(1 / sum(1 / (std.error) ^ 2))) %>% 
  mutate(survey = "YCLS summary",
         psv_type = "retrospective") %>% 
  select(survey, everything())

combined_retrospective_estimates <- 
  bind_rows(psv_retrospective_estimates, abp_retrospective_estimates, 
            ycls_retrospective_estimates) %>% 
  select(survey, Z_psv, outcome_group, estimate) %>% 
  pivot_wider(names_from = survey, values_from = estimate) %>% 
  mutate(abp_diff = `YCLS summary` - `Aronow et al. (2019)`,
         psv_diff = `YCLS summary` - `Press et al. (2013)`,
         abp_ratio = `YCLS summary`/`Aronow et al. (2019)`,
         psv_ratio = `YCLS summary`/`Press et al. (2013)`) %>% 
  gather(survey, estimate, `Press et al. (2013)`:psv_diff) 

combined_retrospective_ses <- 
  bind_rows(psv_retrospective_estimates, abp_retrospective_estimates, 
            ycls_retrospective_estimates) %>% 
  select(survey, Z_psv, outcome_group, std.error) %>% 
  pivot_wider(names_from = survey, values_from = std.error) %>% 
  mutate(abp_diff = sqrt(`Aronow et al. (2019)`^2 + `YCLS summary`^2),
         psv_diff = sqrt(`Press et al. (2013)`^2 + `YCLS summary`^2)) %>% 
  gather(survey, std.error, `Press et al. (2013)`:psv_diff) 

## Table of output 
retrospective_summary <- 
  left_join(combined_retrospective_estimates, combined_retrospective_ses) %>% 
  mutate(p = 2*(1-pnorm(abs(estimate/std.error))), 
         entry = make_entry(est = estimate, se = std.error, p = p, digits = 2),
         outcome_group = ifelse(outcome_group == "Ethical Strike",
                                "Ethical", "Approve")) %>% 
  select(Z_psv, entry, abp_ratio, psv_ratio, outcome_group, survey) %>%
  spread(survey, entry) %>% 
  select(Z_psv, outcome_group, `YCLS summary`, `Press et al. (2013)`, psv_diff, 
         psv_ratio, `Aronow et al. (2019)`, abp_diff, abp_ratio) %>% 
  arrange(outcome_group)

retrospective_summary %>% 
  filter(outcome_group == "Ethical") %>% 
  select(-Z_psv) %>% 
  xtable() %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0("appendix_table_a8b_top.tex"))

retrospective_summary %>% 
  filter(outcome_group == "Approve") %>% 
  select(-Z_psv) %>% 
  xtable() %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0("appendix_table_a8b_btm.tex"))

###-------- Study 9: Attitudes towards immigrants ------------

## Reformat conjoint experiment to long format: 
extract_levels <- function(attributes_vector) {
  Gender <- c("Female", "Male")
  Education <- c(
    "No formal education",
    "Equivalent to completing fourth grade in the U.S.",
    "Equivalent to completing eighth grade in the U.S.",
    "Equivalent to completing high school in the U.S.",
    "Equivalent to completing two years at college in the U.S.",
    "Equivalent to completing a college degree in the U.S.",
    "Equivalent to completing a graduate degree in the U.S."
  )
  
  `Country of Origin` <-
    c(
      "Mexico",
      "Somalia",
      "Iraq",
      "Sudan",
      "China",
      "Philippines",
      "France",
      "Germany",
      "Poland",
      "India"
    )
  
  Profession <-
    c(
      "Research scientist",
      "Child care provider",
      "Nurse",
      "Teacher",
      "Computer programmer",
      "Waiter",
      "Financial analyst",
      "Doctor",
      "Construction worker",
      "Janitor",
      "Gardener"
    )
  
  `Job plans` <- 
    c(
      "Has a contract with a U.S. employer",
      "Does not have a contract with a U.S. employer, but has done job interviews",
      "Will look for work after arriving in the U.S.",
      "Has no plans to look for work at this time"
    )
  
  `Job experience` <- 
    c(
      "No job training or prior experience",
      "One to two years",
      "Three to five years", 
      "More than five years"
    )
  
  `Language Skills` <-
    c(
      "During admission interview, this applicant spoke through an interpreter",
      "During admission interview, this applicant spoke fluent English",
      "During admission interview, this applicant spoke broken English",
      "During admission interview, this applicant tried to speak English but was unable"
    )
  
  `Prior Entry` <-
    c(
      "Never been to the U.S.",
      "Has visited the U.S. many times before on tourist visas",
      "Entered the U.S. once before on a tourist visa",
      "Entered the U.S. once before without legal authorization",
      "Spent six months with family members in the U.S."
    )
  
  `Reason for Application` <-
    c(
      "Escape political/religious persecution",
      "Reunite with family members already in U.S.",
      "Seek better job in U.S."
    )
  
  tibble(
    gender_ = Gender[str_detect(attributes_vector, Gender)],
    education_ = Education[str_detect(attributes_vector, Education)],
    language_ = `Language Skills`[str_detect(attributes_vector, `Language Skills`)],
    country_ = `Country of Origin`[str_detect(attributes_vector, `Country of Origin`)],
    profession_ = Profession[str_detect(attributes_vector, Profession)],
    experience_ = `Job experience`[str_detect(attributes_vector, `Job experience`)], 
    employment_ = `Job plans`[str_detect(attributes_vector, `Job plans`)], 
    reason_ = `Reason for Application`[str_detect(attributes_vector, `Reason for Application`)],
    trips_ = `Prior Entry`[str_detect(attributes_vector, `Prior Entry`)]
  )
  
}


cand_1_1 <-
  ycls_dat %>%
  filter(week == "YCLS Week 8") %>% 
  mutate(conjoint_chosen = as.numeric(conjoint_1_choice == "Immigrant 1"),
         conjoint_rating = conjoint_1_imm_1) %>%
  unite(col = "attributes_vector", starts_with("F-1-1-")) %>%
  group_by(admin_ResponseId, conjoint_chosen, conjoint_rating) %>%
  do(extract_levels(attributes_vector = .$attributes_vector))

cand_1_2 <-
  ycls_dat %>%
  filter(week == "YCLS Week 8") %>% 
  mutate(conjoint_chosen = as.numeric(conjoint_1_choice == "Immigrant 2"),
         conjoint_rating = conjoint_1_imm_2) %>%
  unite(col = "attributes_vector", starts_with("F-1-2-")) %>%
  group_by(admin_ResponseId, conjoint_chosen, conjoint_rating) %>%
  do(extract_levels(attributes_vector = .$attributes_vector))


cand_2_1 <-
  ycls_dat %>%
  filter(week == "YCLS Week 8") %>% 
  mutate(conjoint_chosen = as.numeric(conjoint_2_choice == "Immigrant 1"),
         conjoint_rating = conjoint_2_imm_1) %>%
  unite(col = "attributes_vector", starts_with("F-2-1-")) %>%
  group_by(admin_ResponseId, conjoint_chosen, conjoint_rating) %>%
  do(extract_levels(attributes_vector = .$attributes_vector))

cand_2_2 <-
  ycls_dat %>%
  filter(week == "YCLS Week 8") %>% 
  mutate(conjoint_chosen = as.numeric(conjoint_2_choice == "Immigrant 2"),
         conjoint_rating = conjoint_2_imm_2) %>%
  unite(col = "attributes_vector", starts_with("F-2-2-")) %>%
  group_by(admin_ResponseId, conjoint_chosen, conjoint_rating) %>%
  do(extract_levels(attributes_vector = .$attributes_vector))


cand_3_1 <-
  ycls_dat %>%
  filter(week == "YCLS Week 8") %>% 
  mutate(conjoint_chosen = as.numeric(conjoint_3_choice == "Immigrant 1"),
         conjoint_rating = conjoint_3_imm_1) %>%
  unite(col = "attributes_vector", starts_with("F-3-1-")) %>%
  group_by(admin_ResponseId, conjoint_chosen, conjoint_rating) %>%
  do(extract_levels(attributes_vector = .$attributes_vector))

cand_3_2 <-
  ycls_dat %>%
  filter(week == "YCLS Week 8") %>% 
  mutate(conjoint_chosen = as.numeric(conjoint_3_choice == "Immigrant 2"),
         conjoint_rating = conjoint_3_imm_2) %>%
  unite(col = "attributes_vector", starts_with("F-3-2-")) %>%
  group_by(admin_ResponseId, conjoint_chosen, conjoint_rating) %>%
  do(extract_levels(attributes_vector = .$attributes_vector))

cand_4_1 <-
  ycls_dat %>%
  filter(week == "YCLS Week 8") %>% 
  mutate(conjoint_chosen = as.numeric(conjoint_4_choice == "Immigrant 1"),
         conjoint_rating = conjoint_4_imm_1) %>%
  unite(col = "attributes_vector", starts_with("F-4-1-")) %>%
  group_by(admin_ResponseId, conjoint_chosen, conjoint_rating) %>%
  do(extract_levels(attributes_vector = .$attributes_vector))

cand_4_2 <-
  ycls_dat %>%
  filter(week == "YCLS Week 8") %>% 
  mutate(conjoint_chosen = as.numeric(conjoint_4_choice == "Immigrant 2"),
         conjoint_rating = conjoint_4_imm_2) %>%
  unite(col = "attributes_vector", starts_with("F-4-2-")) %>%
  group_by(admin_ResponseId, conjoint_chosen, conjoint_rating) %>%
  do(extract_levels(attributes_vector = .$attributes_vector))

cand_5_1 <-
  ycls_dat %>%
  filter(week == "YCLS Week 8") %>% 
  mutate(conjoint_chosen = as.numeric(conjoint_5_choice == "Immigrant 1"),
         conjoint_rating = conjoint_5_imm_1) %>%
  unite(col = "attributes_vector", starts_with("F-5-1-")) %>%
  group_by(admin_ResponseId, conjoint_chosen, conjoint_rating) %>%
  do(extract_levels(attributes_vector = .$attributes_vector))

cand_5_2 <-
  ycls_dat %>%
  filter(week == "YCLS Week 8") %>% 
  mutate(conjoint_chosen = as.numeric(conjoint_5_choice == "Immigrant 2"),
         conjoint_rating = conjoint_5_imm_2) %>%
  unite(col = "attributes_vector", starts_with("F-5-2-")) %>%
  group_by(admin_ResponseId, conjoint_chosen, conjoint_rating) %>%
  do(extract_levels(attributes_vector = .$attributes_vector))

conjoint_replication <-
  bind_rows(
    cand_1_1 = cand_1_1,
    cand_1_2 = cand_1_2,
    cand_2_1 = cand_2_1,
    cand_2_2 = cand_2_2,
    cand_3_1 = cand_3_1,
    cand_3_2 = cand_3_2,
    cand_4_1 = cand_4_1,
    cand_4_2 = cand_4_2,
    cand_5_1 = cand_5_1,
    cand_5_2 = cand_5_2,
    .id = "candidate"
  ) %>%
  ungroup %>%
  left_join(ycls_dat %>%
              filter(week == "YCLS Week 8")) %>%
  select(
    starts_with("admin_"),
    starts_with("x_"),
    starts_with("conjoint_"),
    ends_with("_")
  )


## Original conjoint experiment
conjoint_original <-
  read_dta("conjoint_original.dta") %>%
  mutate(
    education_ = as_factor(FeatEd),
    gender_ = as_factor(FeatGender),
    country_ = as_factor(FeatCountry),
    reason_ = as_factor(FeatReason),
    profession_ = as_factor(FeatJob),
    experience_ = as_factor(FeatExp),
    employment_ = as_factor(FeatPlans),
    trips_ = as_factor(FeatTrips),
    language_ = as_factor(FeatLang))

## Relabel replication for consistency w/ original
conjoint_replication <-
  conjoint_replication %>%
  mutate(
    education_ =
      case_when(
        education_ == "Equivalent to completing two years at college in the U.S." ~
          "Equivalent to completing two years of college in the U.S.",
        TRUE ~ education_
      ),
    education_ = gsub(
      pattern = "\\.",
      replacement = "",
      x = education_
    ),
    gender_ = tolower(gender_),
    reason_ =
      case_when(
        reason_ == "Reunite with family members already in U.S." ~
          "Reunite with family members already in the U.S.",
        TRUE ~ reason_
      ),
    employment_ =
      case_when(
        employment_ == "Does not have a contract with a U.S. employer, but has done job interviews" ~
          "Does not have a contract with a U.S. employer but has done job interviews",
        TRUE ~ employment_
      ),
    trips_ =
      case_when(
        trips_ == "Entered the U.S. once before on a tourist visa" ~
          "Entered U.S. once before on a tourist visa",
        trips_ == "Spent six months with family members in the U.S." ~
          "Spent six months with family members in the U.S",
        TRUE ~ trips_
      ),
    language_ =
      case_when(
        language_ == "During admission interview, this applicant spoke fluent English" ~
          "During the admission interview, this applicant spoke fluent English",
        language_ == "During admission interview, this applicant spoke through an interpreter" ~
          "During admission interview, this applicant spoke [language] and used an interpreter",
        TRUE ~ language_
      ),
    experience_ =
      case_when(
        experience_ == "One to two years" ~ "One or two years of job training and experience",
        experience_ == "Three to five years" ~ "Three to five years of job training and experience",
        experience_ == "More than five years" ~ "More than five years of job training and experience",
        TRUE ~ experience_
      ),
    education_ = factor(
      education_,
      levels = levels(conjoint_original$education_),
      labels = c(
        "No formal",
        "4th grade",
        "8th grade",
        "High school",
        "Two-year college",
        "College degree",
        "Graduate degree"
      )
    ),
    gender_ = factor(
      gender_,
      levels = levels(conjoint_original$gender_),
      labels = c("Female", "Male")
    ),
    country_ = factor(country_, levels = levels(conjoint_original$country_)),
    reason_ = factor(
      reason_,
      levels = levels(conjoint_original$reason_),
      labels = c("Reunite with family",
                 "Seek better job",
                 "Escape persecution")
    ),
    profession_ = factor(profession_, levels = levels(conjoint_original$profession_)),
    experience_ = factor(
      experience_,
      levels = levels(conjoint_original$experience_),
      labels = c("None", "1-2 years", "3-5 years",
                 "5+ years")
    ),
    employment_ = factor(
      employment_,
      levels = levels(conjoint_original$employment_),
      labels = c(
        "Contract with employer",
        "Interviews with employer",
        "Will look for work",
        "No plans to look for work"
      )
    ),
    trips_ = factor(
      trips_,
      levels = levels(conjoint_original$trips_),
      labels = c(
        "Never",
        "Once as tourist",
        "Many times as tourist",
        "Six months with family",
        "Once w/o authorization"
      )
    ),
    language_ = factor(
      language_,
      levels = levels(conjoint_original$language_),
      labels = c(
        "Fluent English",
        "Broken English",
        "Tried English but unable",
        "Used interpreter"
      )
    )
  )

conjoint_original <-
  conjoint_original %>%
  mutate(
    education_ = factor(
      education_,
      labels = c(
        "No formal",
        "4th grade",
        "8th grade",
        "High school",
        "Two-year college",
        "College degree",
        "Graduate degree"
      )
    ),
    gender_ = factor(gender_, labels = c("Female", "Male")),
    
    reason_ = factor(
      reason_,
      labels = c("Reunite with family",
                 "Seek better job",
                 "Escape persecution")
    ),
    experience_ = factor(
      experience_,
      labels = c("None", "1-2 years", "3-5 years",
                 "5+ years")
    ),
    employment_ = factor(
      employment_,
      labels = c(
        "Contract with employer",
        "Interviews with employer",
        "Will look for work",
        "No plans to look for work"
      )
    ),
    trips_ = factor(
      trips_,
      labels = c(
        "Never",
        "Once as tourist",
        "Many times as tourist",
        "Six months with family",
        "Once w/o authorization"
      )
    ),
    language_ = factor(
      language_,
      labels = c(
        "Fluent English",
        "Broken English",
        "Tried English but unable",
        "Used interpreter"
      )
    )
  )


## Estimate AMCEs and combine
original_ests <-
  lm_robust(
    Chosen_Immigrant ~ gender_ + education_ + language_ + country_ + 
      profession_ + experience_ + employment_ + reason_ + trips_,
    clusters = CaseID,
    weights = weight2,
    data = conjoint_original
  ) %>%
  tidy


replication_ests <-
  lm_robust(
    conjoint_chosen ~ gender_ + education_ + language_ + country_ + 
      profession_ + experience_ + employment_ + reason_ + trips_,
    clusters = admin_ResponseId,
    data = conjoint_replication
  ) %>%
  tidy

est_study_9 <- 
  bind_rows(`Hainmueller & Hopkins (2015)` = original_ests,
          `Week 8` = replication_ests,
          .id = "survey") %>%
  filter(term != "(Intercept)") %>%
  separate(term, into = c("attribute", "level"), sep = "_") %>%
  mutate(
    attribute_ref =
      case_when(
        attribute == "gender" ~ "Gender (reference: Female applicant)",
        attribute == "education" ~ "Education (reference: No formal education)",
        attribute == "language" ~ "Language (reference: Applicant spoke fluent English)",
        attribute == "country" ~ "Country of origin (reference: Germany)",
        attribute == "profession" ~ "Profession (reference: Janitor)",
        attribute == "reason" ~ "Application reason (reference: Reunite with family members already in the U.S.)",
        attribute == "trips" ~ "Prior trips to U.S. (reference: Never been to the U.S.)",
        attribute == "experience" ~ "Job experience (reference: No job training or prior experience)",
        attribute == "employment" ~ "Jobs plans (reference: Contract with a U.S. employer)"
      ),
    attribute_ref = factor(
      attribute_ref,
      levels = c(
        "Gender (reference: Female applicant)",
        "Education (reference: No formal education)",
        "Language (reference: Applicant spoke fluent English)",
        "Country of origin (reference: Germany)",
        "Profession (reference: Janitor)",
        "Job experience (reference: No job training or prior experience)",
        "Jobs plans (reference: Contract with a U.S. employer)",
        "Application reason (reference: Reunite with family members already in the U.S.)",
        "Prior trips to U.S. (reference: Never been to the U.S.)"
      )
    ),
    level = factor(level,
                   levels = rev(
                     c(
                       levels(conjoint_original$gender_)[-1],
                       levels(conjoint_original$education_)[-1],
                       levels(conjoint_original$language_)[-1],
                       levels(conjoint_original$country_)[-1],
                       levels(conjoint_original$profession_)[-1],
                       levels(conjoint_original$experience_)[-1],
                       levels(conjoint_original$employment_)[-1],
                       levels(conjoint_original$reason_)[-1],
                       levels(conjoint_original$trips_)[-1]
                     )
                   )),
    survey = factor(survey, levels = rev(c("Hainmueller & Hopkins (2015)", 
                                           "Week 8"))),
    study = ifelse(survey == "Week 8", "Replication", "Pre-COVID")
  )


### ----------- Export Figure A.9 for Online Appendix ------
g <-
  ggplot(est_study_9,
         aes(
           estimate,
           level,
           group = survey,
           color = survey,
           shape = survey
         )) +
  geom_vline(xintercept = 0,
             lty = 2,
             size = 0.5) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .9),
    size = 1
  ) +
  geom_point(position = position_dodgev(height = .9), size = 2) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_color_grey("", start = 0.2, end = 0.6) +
  scale_shape_manual("", values = c(15, 16)) +
  xlab("Average treatment effect estimate") +
  ylab("") +
  ggtitle("") +
  theme_ycls() + 
  theme(
    strip.text = element_text(size = 10.28),
    strip.text.x = element_text(size = 10.28),
    plot.margin = unit(c(0,0,0,0), "lines")
  ) + 
  geom_text(
    data = (. %>% filter(level ==  "Male")),
    aes(label = study, x = conf.high + 0.1), # c(0.09, 0.06)), 
    position = position_dodgev(height = 1.5),
    size = 3, family = "serif") 

g

g %>%
  ggsave(
    filename = "appendix_fig_a9.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * (nrow(g$data)*.4)
  )

###-------- Study 10: Fake news corrections ------------

## Combine replication dataset w/ original 
pwk_original <-
  read_rds("pwk_original.rds") %>%
  mutate(
    treatment_2 = case_when(
      treatment %in% c("green correction", "wapost correction") ~ "correction",
      TRUE ~ as.character(treatment)
    ),
    treatment_2 = factor(treatment_2, levels = c("no correction", "correction")),
    fakenews_scaramucci_Y_std = glass_delta(
      Y = 6 - fakenews_scaramucci_Y,
      Z = treatment_2,
      level = "no correction"
    ),
    fakenews_pizza_Y_std = glass_delta(
      Y = 6 - fakenews_pizza_Y,
      Z = treatment_2,
      level = "no correction"
    ),
    fakenews_podesta_Y_std = glass_delta(
      Y = 6 - fakenews_podesta_Y,
      Z = treatment_2,
      level = "no correction"
    ),
    fakenews_pedo_Y_std = glass_delta(
      Y = 6 - fakenews_pedo_Y,
      Z = treatment_2,
      level = "no correction"
    ),
    fakenews_obama_Y_std = glass_delta(
      Y = 6 - fakenews_obama_Y,
      Z = treatment_2,
      level = "no correction"
    ),
    fakenews_vermont_Y_std = glass_delta(
      Y = 6 - fakenews_vermont_Y,
      Z = treatment_2,
      level = "no correction"
    )
  )

pwk_replication <-
  ycls_dat %>%
  filter(week == "YCLS Week 4") %>% 
  mutate(
    fakenews_vermont_Z2 = if_else(
      fakenews_vermont_Z == "no correction",
      "no correction",
      "correction"
    ),
    fakenews_vermont_Z2 = factor(fakenews_vermont_Z2, levels = c("no correction", "correction")),
    fakenews_scaramucci_Y_std = glass_delta(
      Y = 6 - fakenews_scaramucci_Y,
      Z =  fakenews_scaramucci_Z,
      level = "no correction"
    ),
    fakenews_pizza_Y_std = glass_delta(
      Y = 6 - fakenews_pizza_Y,
      Z = fakenews_pizza_Z,
      level = "no correction"
    ),
    fakenews_podesta_Y_std = glass_delta(
      Y = 6 - fakenews_podesta_Y,
      Z =  fakenews_podesta_Z,
      level = "no correction"
    ),
    fakenews_pedo_Y_std = glass_delta(
      Y = 6 - fakenews_pedo_Y,
      Z =  fakenews_pedo_Z,
      level = "no correction"
    ),
    fakenews_obama_Y_std = glass_delta(
      Y = 6 - fakenews_obama_Y,
      Z =  fakenews_obama_Z,
      level = "no correction"
    ),
    fakenews_vermont_Y_std = glass_delta(
      Y = 6 - fakenews_vermont_Y,
      Z =  fakenews_vermont_Z2,
      level = "no correction"
    )
  )


## Replication estimates:
fit_1 <-
  lm_robust(fakenews_scaramucci_Y_std ~ fakenews_scaramucci_Z,
            data = pwk_replication)

fit_2 <-
  lm_robust(fakenews_pizza_Y_std ~ fakenews_pizza_Z,
            data = pwk_replication)

fit_3 <-
  lm_robust(fakenews_podesta_Y_std ~ fakenews_podesta_Z,
            data = pwk_replication)

fit_4 <-
  lm_robust(fakenews_pedo_Y_std ~ fakenews_pedo_Z,
            data = pwk_replication)

fit_5 <-
  lm_robust(fakenews_obama_Y_std ~ fakenews_obama_Z,
            data = pwk_replication)

fit_6 <-
  lm_robust(fakenews_vermont_Y_std ~ fakenews_vermont_Z2,
            data = pwk_replication)

replication_est <-
  list(fit_1, fit_2, fit_3, fit_4, fit_5, fit_6) %>%
  map_df(tidy) %>%
  filter(term != "(Intercept)")

## Original estimates: 
fit_7 <-
  lm_robust(fakenews_scaramucci_Y_std ~ treatment_2, data = pwk_original)

fit_8 <-
  lm_robust(fakenews_podesta_Y_std ~ treatment_2, data = pwk_original)

fit_9 <-
  lm_robust(fakenews_pizza_Y_std ~ treatment_2, data = pwk_original)

fit_10 <-
  lm_robust(fakenews_pedo_Y_std ~ treatment_2, data = pwk_original)

fit_11 <-
  lm_robust(fakenews_obama_Y_std ~ treatment_2, data = pwk_original)

fit_12 <-
  lm_robust(fakenews_vermont_Y_std ~ treatment_2, data = pwk_original)

original_est <-
  list(fit_7, fit_8, fit_9, fit_10, fit_11, fit_12) %>%
  map_df(tidy) %>%
  filter(term != "(Intercept)")

## Combine and prep for plot
est_study_10 <-
  bind_rows(`Porter et al. (2018)` = original_est,
            `Week 4` = replication_est,
            .id = "survey") %>%
  mutate(
    outcome_label =
      factor(
        outcome,
        levels = c(
          "fakenews_pizza_Y_std",
          "fakenews_podesta_Y_std",
          "fakenews_vermont_Y_std",
          "fakenews_scaramucci_Y_std",
          "fakenews_pedo_Y_std",
          "fakenews_obama_Y_std"
        ),
        labels = c(
          "DC pizza restaurant concealed sex\ndungeon used by Democratic elites",
          "John Podesta implicated in\ndisappearance of Madeleine McCann",
          "Russian government hacks\nVermont power plant",
          "Anthony Scaramucci subject of\nSenate Russia investigation",
          "President Trump orders\ncrackdown on sex trafficking",
          "President Obama fakes\nbirth certificate"
        )
      ),
    survey = factor(survey, levels = c("Week 4", "Porter et al. (2018)")),
    study = factor(
      survey,
      levels = c("Week 4", "Porter et al. (2018)"),
      labels = c("Replication", "Pre-COVID")
    ),
    estimate_type = "Standardized"
  )

### ----------- Export Figure A.10 for Online Appendix ------
g <-
  ggplot(est_study_10,
         aes(
           estimate,
           outcome_label,
           group = survey,
           color = survey,
           shape = survey
         )) +
  geom_vline(xintercept = 0,
             lty = 2,
             size = 0.5) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.5),
    size = 1
  ) +
  geom_point(position = position_dodgev(height = 0.5), size = 2.5) +
  scale_color_grey("", start = 0.2, end = 0.6) +
  scale_shape_manual("", values = c(15, 16)) +
  xlab("Average treatment effect estimate (standard units)") +
  ylab("") +
  ggtitle("") +
  theme_ycls() +
  geom_text(
    data = (
      . %>% filter(outcome_label == "President Obama fakes\nbirth certificate")
    ),
    aes(x = c(0.6, 0.6), label = study),
    size = 3,
    position = position_dodgev(height = 0.5)
  )

g

g %>%
  ggsave(
    filename = "appendix_fig_a10.pdf",
    .,
    width = 6.5,
    height = fixed + spacing * nrow(g$data)
  )

###-------- Study 11: Inequality and system justification ------------

tw_original <-
  read_rds("tw_original.rds") %>% 
  mutate(
    tw_scale_1_s = glass_delta(Y = tw_scale_1, Z = tw_Z,
                               level = "control"),
    tw_scale_2_s = glass_delta(Y = tw_scale_2, Z = tw_Z,
                               level = "control"),
    tw_scale_3_s = glass_delta(Y = tw_scale_3, Z = tw_Z,
                               level = "control"),
    tw_manip_1_s = glass_delta(Y = tw_manip_1, Z = tw_Z,
                               level = "control"),
    tw_manip_2_s = glass_delta(Y = tw_manip_2, Z = tw_Z,
                               level = "control")
  )

tw_replication <-
  ycls_dat %>%
  filter(week == "YCLS Week 2") %>% 
  mutate(
    tw_scale_1_s = glass_delta(Y = tw_scale_1, Z = tw_Z,
                               level = "control"),
    tw_scale_2_s = glass_delta(Y = tw_scale_2, Z = tw_Z,
                               level = "control"),
    tw_scale_3_s = glass_delta(Y = tw_scale_3, Z = tw_Z,
                               level = "control"),
    tw_manip_1_s = glass_delta(Y = tw_manip_1, Z = tw_Z,
                               level = "control"),
    tw_manip_2_s = glass_delta(Y = tw_manip_2, Z = tw_Z,
                               level = "control"),
    tw_Z = relevel(tw_Z, ref = "control")
  ) %>% 
  filter(tw_Z != "pure_control")

## Estimate treatment effects and combine: 
fit_1a <-
  lm_robust(tw_manip_1 ~ tw_Z + x_pid7, data = tw_replication)

fit_2a <-
  lm_robust(tw_manip_2 ~ tw_Z + x_pid7, data = tw_replication)

fit_1b <-
  lm_robust(tw_manip_1_s ~ tw_Z + x_pid7, data = tw_replication)

fit_2b <-
  lm_robust(tw_manip_2_s ~ tw_Z + x_pid7, data = tw_replication)

fit_3 <-
  lm_robust(tw_scale_1_s ~ tw_Z + x_pid7, data = tw_replication)

fit_4 <-
  lm_robust(tw_scale_2_s ~ tw_Z + x_pid7, data = tw_replication)

fit_5 <-
  lm_robust(tw_scale_3_s ~ tw_Z + x_pid7, data = tw_replication)

fit_6a <-
  lm_robust(tw_manip_1 ~ tw_Z + pid_7n, data = tw_original)

fit_7a <-
  lm_robust(tw_manip_2 ~ tw_Z + pid_7n, data = tw_original)

fit_6b <-
  lm_robust(tw_manip_1_s ~ tw_Z + pid_7n, data = tw_original)

fit_7b <-
  lm_robust(tw_manip_2_s ~ tw_Z + pid_7n, data = tw_original)

fit_8 <-
  lm_robust(tw_scale_1_s ~ tw_Z + pid_7n, data = tw_original)

fit_9 <-
  lm_robust(tw_scale_2_s ~ tw_Z + pid_7n, data = tw_original)

fit_10 <-
  lm_robust(tw_scale_3_s ~ tw_Z + pid_7n, data = tw_original)

ycls_est <-
  list(fit_1a,
       fit_2a,
       fit_1b,
       fit_2b,
       fit_3,
       fit_4,
       fit_5) %>%
  map_df(tidy)

original_est <-
  list(fit_6a,
       fit_7a,
       fit_6b,
       fit_7b,
       fit_8,
       fit_9,
       fit_10) %>%
  map_df(tidy)

## Combine estimates and prep for plotting
est_study_11 <-
  bind_rows(`Week 2` = ycls_est,
            `Trump & White (2018)` = original_est,
            .id = "survey") %>%
  filter(!term %in% c("(Intercept)", "x_pid7", "pid_7n")) %>%
  mutate(
    survey = factor(survey, levels = c("Week 2", "Trump & White (2018)")),
    study = ifelse(survey == "Week 2", "Replication", "Pre-COVID"),
    outcome_group =
      case_when(
        outcome %in% c("tw_scale_1_s", "tw_scale_2_s", "tw_scale_3_s") ~ "Outcome Scales (standard units)",
        outcome %in% c("tw_manip_1", "tw_manip_2") ~ "Manipulation Checks (binary)",
        outcome %in% c("tw_manip_1_s", "tw_manip_2_s") ~ "Manipulation Checks (standard units)"
      ),
    outcome_factor =
      factor(
        outcome,
        levels = c(
          "tw_manip_1",
          "tw_manip_2",
          "tw_manip_1_s",
          "tw_manip_2_s",
          "tw_scale_1_s",
          "tw_scale_2_s",
          "tw_scale_3_s"
        ),
        labels = c(
          "Inequality increased",
          "Richs' share of income changed",
          "Inequality increased",
          "Richs' share of income changed",
          "System Justification",
          "Institutional Trust",
          "Econ. System Justification"
        )
      ),
    estimate_type = case_when(
     str_detect(outcome, "_s") ~ "Standardized", 
     TRUE ~ "Unstandardized"
    )
  )

### ----------- Export Figure A.11 for Online Appendix ------
g <-
  est_study_11 %>%
  filter(
    outcome_group %in% c(
      "Manipulation Checks (binary)",
      "Outcome Scales (standard units)"
    )
  ) %>%
  ggplot(.,
         aes(
           estimate,
           outcome_factor,
           group = survey,
           color = survey,
           shape = survey
         )) +
  geom_vline(xintercept = 0,
             lty = 2,
             size = 0.5) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.5),
    size = 1
  ) +
  geom_point(position = position_dodgev(height = 0.5), size = 2.5) +
  facet_rep_wrap( ~ outcome_group, ncol = 1, scales = "free") +
  scale_color_grey("", start = 0.2, end = 0.6) +
  scale_shape_manual("", values = c(15, 16)) +
  xlab("Average treatment effect estimate") +
  ylab("") +
  theme_ycls() +
  theme(
    axis.line.x = element_line(size = 0.5),
    panel.spacing.x = unit(-5, "lines"),
    strip.text.x = element_text(hjust = 0.5)
  ) +
  coord_capped_cart(bottom = "none") +
  geom_text(
    data = (. %>% filter(
      outcome_factor ==  "Richs' share of income changed"
    )),
    aes(label = study, x = c(0.3, 0.3)),
    position = position_dodgev(height = 0.5),
    size = 3,
    family = "serif"
  )


g

g %>%
  ggsave(
    filename = "appendix_fig_a11.pdf",
    .,
    width = 6.5,
    height = fixed + 1 + spacing * nrow(g$data)
  )


###-------- Study 12: Trust in government and redistribution ------------
peyton_original <-
  read_csv("peyton_original.csv") %>%
  mutate(
    Z_trust = factor(
      Z,
      levels = c("Corrupt", "Placebo", "Honest"),
      labels = c("Corrupt", "Control", "Honest")
    ),
    Z_trust_n = case_when(
      Z_trust == "Corrupt" ~ 0,
      Z_trust == "Control" ~ 0.5,
      Z_trust == "Honest" ~ 1
    ),
    D = glass_delta(Y = fmptrust_index, Z = Z_trust,
                    level = "Control"),
    Y = glass_delta(Y = fedspend_index, Z = Z_trust,
                    level = "Control")
  ) %>%
  filter(study_type == "Politics") 

peyton_rep <-
  ycls_dat %>% 
  filter(str_detect(week, "Week 9")) %>%
  mutate(
    Z_trust = factor(
      peyton_Z,
      levels = c("corrupt", "placebo",
                 "honest"),
      labels = c("Corrupt", "Control", "Honest")
    ),
    Z_trust_n = case_when(
      Z_trust == "Corrupt" ~ 0,
      Z_trust == "Control" ~ 0.5,
      Z_trust == "Honest" ~ 1
    ),
    D = glass_delta(Y = peyton_trust_index, Z = Z_trust,
                    level = "Control"),
    Y = glass_delta(Y = peyton_fedspend_index, Z = Z_trust,
                    level = "Control")
  )

combined_dat <-
  bind_rows(
    `MTurk sample, Jun 2014` = peyton_original %>% filter(experiment == "MTurk_Pols14"),
    `Qualtrics sample, Sep 2014` = peyton_original %>% filter(experiment == "GenPop_Pols14"),
    `MTurk sample, Mar 2017` = peyton_original %>% filter(experiment == "MTurk_Pols17"),
    `Week 9 Replication` = peyton_rep,
    .id = "survey"
  ) %>%
  filter(!is.na(Z_trust))

## Estimate treatment effects and combine: 

# First stage fx of Z on D
trustgov <-
  combined_dat %>%
  group_by(survey) %>%
  do(tidy(lm_robust(
    D ~ Z_trust_n, 
    data = .
  ))) 

# Reduced form fx of Z on Y
fedspend <-
  combined_dat %>% 
  group_by(survey) %>%
  do(tidy(lm_robust(
    Y ~ Z_trust_n, 
    data = .
  ))) 

est_study_12 <-
  bind_rows(trustgov, fedspend) %>%
  filter(term != "(Intercept)") %>% 
  mutate(
    outcome = factor(
      outcome,
      levels = c("D", "Y"),
      labels = c("Trust in Government",
                 "Support for Redistribution")
    ),
    survey = factor(survey, levels = rev(
      c(
        "MTurk sample, Jun 2014",
        "Qualtrics sample, Sep 2014",
        "MTurk sample, Mar 2017",
        "Week 9 Replication"
      )
    )),
    study = ifelse(survey == "Week 9 Replication", "Replication", "Pre-COVID"),
    estimate_type = "Standardized"
  )


### ----------- Export Figure A.12 for Online Appendix ------
g <-
  est_study_12 %>%
  ggplot(.,
         aes(
           x = estimate,
           y = survey,
           color = study,
           shape = study
         )) +
  geom_vline(
    xintercept = 0,
    lty = 2,
    size = 0.5
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.5), size = 1
  ) +
  geom_point(position = position_dodgev(height = 0.5), size = 2.5 
  ) +
  facet_rep_wrap(~ outcome, ncol = 1) + 
  scale_color_grey("", start = 0.2, end = 0.6) +
  scale_shape_manual("", values = c(15, 16)) +
  xlab("Average treatment effect estimate (standard units)") +
  ylab("") +
  theme_ycls() + 
  theme(axis.line.x = element_line(size = 0.5),
        strip.text.x = element_text(hjust = 0.5)) + 
  coord_capped_cart(bottom = "none") 

g


g %>%
  ggsave(
    filename = "appendix_fig_a12.pdf",
    .,
    width = 6.5,
    height = fixed + 1 + spacing * nrow(g$data)
  )



###-------- Combine and Export Summary Effect Size Estimates --------

conjoint_levs <- unique(est_study_9$level)

combined_dat <-
  bind_rows(
    `Hyman & Sheatsley (1950)` = est_study_1,
    `Tversky & Kaheneman (1981) - Cheap/Expensive` = est_study_2,
    `Tversky & Kaheneman (1981) - Gain/Loss` = est_study_3,
    `Smith (1987)` = est_study_4 %>%
      filter(type == "Experimental", !(survey %in% paste("GSS", 1987:2018))),
    `Druckman (2001)` = est_study_5,
    `Gilens (2001)` = est_study_6,
    `Knobe (2003)` = est_study_7,
    `Press et al. (2013)` = bind_rows(est_study_8a, est_study_8b),
    `Hainmueller & Hopkins (2015)` = est_study_9,
    `Porter et al. (2018)` = est_study_10,
    `Trump & White (2018)` = est_study_11,
    `Peyton (2020)` = est_study_12,
    .id = "study_group"
  ) %>%
  mutate(
    study_group_detail = case_when(
      outcome == "psv_prefer_binary_s" &
        term == "90/45" ~ "Prefer Nukes (90/45)",
      outcome == "psv_approve_nuke_binary_s" &
        term == "90/45"  ~ "Approve Nukes (90/45)",
      outcome == "psv_prefer_binary_s" &
        term == "90/70" ~ "Prefer Nukes (90/70)",
      outcome == "psv_approve_nuke_binary_s" &
        term == "90/70"  ~ "Approve Nukes (90/70)",
      outcome == "psv_retrospective_ethical_binary_s" ~ "Ethical Strike (Nuclear)",
      outcome == "psv_retrospective_approve_binary_s" ~ "Approve Strike (Nuclear)",
      AD_Z_party_sure_thing == "Republicans'\n Program" ~ "Republicans' Program",
      AD_Z_party_sure_thing == "Democrats'\n Program" ~ "Democrats' Program",
      AD_Z_party_sure_thing == "Program A" ~ "Program A",
      outcome == "Trust in Government" ~ "Trust in Government",
      outcome == "Support for Redistribution" ~ "Redistribution",
      outcome == "tw_manip_2_s" ~ "Rich's share of income changed",
      outcome == "tw_manip_1_s" ~ "Inequality increased",
      outcome == "tw_scale_1_s" ~ "System Justification",
      outcome == "tw_scale_2_s" ~ "Institutional Trust",
      outcome == "tw_scale_3_s" ~ "Economic System Justification",
      outcome == "fakenews_scaramucci_Y_std" ~ "Scaramucci",
      outcome == "fakenews_podesta_Y_std" ~ "Podesta",
      outcome == "fakenews_pedo_Y_std" ~ "Sex trafficking",
      outcome == "fakenews_vermont_Y_std" ~ "Russian hackers",
      outcome == "fakenews_obama_Y_std" ~ "Obama Birtherism",
      outcome == "fakenews_pizza_Y_std" ~ "D.C. Pizzagate",
      term == "framing_ZCheap" ~ "Cheap v. Expensive Framing",
      term == "AD_Z_loss" ~ "Gain v. Loss Framing",
      study_group == "Smith (1987)" ~ "Welfare v. Aid to Poor",
      study_group == "Knobe (2003)" ~ "Perceived Intentionality",
      study_group == "Hyman & Sheatsley (1950)" ~ "Russian reporters",
      study_group == "Gilens (2001)" ~ "Foreign aid misperceptions",
      !is.na(level) ~ paste(level),
      TRUE ~ paste(study_group)
    ),
    study_group_detail = factor(
      study_group_detail,
      levels = c(
        "Russian reporters",
        "Cheap v. Expensive Framing",
        "Gain v. Loss Framing",
        "Welfare v. Aid to Poor",
        "Republicans' Program",
        "Democrats' Program",
        "Program A",
        "Foreign aid misperceptions",
        "Perceived Intentionality",
        "Approve Nukes (90/45)",
        "Approve Nukes (90/70)",
        "Prefer Nukes (90/45)",
        "Prefer Nukes (90/70)",
        "Ethical Strike (Nuclear)",
        "Approve Strike (Nuclear)",
        "D.C. Pizzagate",
        "Podesta",
        "Russian hackers",
        "Scaramucci",
        "Sex trafficking",
        "Obama Birtherism",
        "Economic System Justification",
        "Institutional Trust",
        "System Justification",
        "Inequality increased",
        "Rich's share of income changed",
        "Redistribution",
        "Trust in Government",
        paste(conjoint_levs)
      )
    )
  ) %>% 
  filter(estimate_type == "Standardized" |
           is.na(estimate_type))


## Export dataset of summary effect size estimates
summary_dat <-
  combined_dat %>%
  ungroup() %>%
  ## Rescale estimates for plotting so DVs w/ negative sign in pre-COVID
  ## studies are presented in terms of absolute value
  mutate(
    estimate = case_when(
      str_detect(outcome, "psv_retrospective") ~ -1 * estimate,
      str_detect(study_group, "Hainmueller|Trump") ~ ifelse(estimate < 0,
                                                            estimate * (-1),
                                                            estimate),
      TRUE ~ estimate
    ),
    conf.low = case_when(
      str_detect(outcome, "psv_retrospective") ~ -1 * conf.high,
      str_detect(study_group, "Hainmueller|Trump") ~ ifelse(estimate < 0,
                                                            conf.high * (-1),
                                                            conf.low),
      TRUE ~ conf.low
    ),
    conf.high = case_when(
      str_detect(outcome, "psv_retrospective") ~ -1 * conf.low,
      str_detect(study_group, "Hainmueller|Trump") ~ ifelse(estimate < 0,
                                                            conf.low * (-1),
                                                            conf.high),
      TRUE ~ conf.high
    )
  ) %>%
  group_by(study_group_detail, study) %>%
  summarise(
    estimate = weighted.mean(estimate, w = 1 / (std.error) ^ 2),
    std.error = sqrt(1 / sum(1 / (std.error) ^ 2)),
    study_group = paste(study_group)
  ) %>%
  mutate(
    conf.low = estimate - 1.96 * std.error,
    conf.high = estimate + 1.96 * std.error,
    statistic = estimate / std.error,
    p.value = 2 * (1 - pnorm(abs(statistic))),
    study = ifelse(study == "Replication", "YCLS summary",
                   "Pre-COVID summary")
  ) %>%
  distinct()

write_rds(summary_dat, "phc_summary.rds")




